1 Descripción del ejercicio

El presente ejercicio se engloba dentro de la materia de Sistemas de Información Geográfica y Ecología Espacial: Aplicaciones del Máster Geoforest de la Universidad de Córdoba.

El trabajo consiste en aplicar las herramientas aprendidas en clase, con lo que los productos finales esperados son:

1.1 Datos de partida

Para la realización del ejercicio partimos de los siguientes datos:

  • limite_yunquera.shp: límite de la zona de estudio (EPSG: 25830).

  • vegetacion_yunquera.shp: mapa de vegetación de la zona de estudio.

  • inventario_pinsapo.shp: puntos de inventario de parcelas de pinsapo (Abies pinsapo).

  • PNOA_MDT05_ETRS89_HU30_1051_LID.tif: modelo digital del terreno con tamaño de píxel de 5x5 metros.

Estos datos se encuentran en la carpeta Cartografia/CartografiaOriginal.

El resto de datos de partida se obtendrán a lo largo de la resolución del ejercicio.

1.2 Zona de estudio

Antes de comenzar a resolver el ejercicio conviene realizar una breve descripción de la zona de estudio.

Yunquera es un municipio perteneciente a la provincia de Málaga (Andalucía), que se sitúa al oeste de la provincia.

Al oeste del municipio de Yunquera, se encuentra el monte Pinar de Yunquera. En este monte encontramos especies arbóreas como Abies pinsapo, Castanea sativa, Cedrus atlantica, Quercus faginea, Quercus ilex subsp. ballota, Pinus halepensis, Pinus nigra y Olea europaea var. europaea. Estos se encuentran formando distintas formaciones entre las que destacan los pinsapares en masas monoespecíficas y pluriespecíficas.

En el siguiente webmap se pueden explorar las diferentes formaciones vegetales presentes en el área de estudio (hacer click en las parcelas para explorar).

yunquera <- read_sf(here('Cartografia/CartografiaOriginal/limite_yunquera.shp'),
                    crs = 25830)
vegetacion <- read_sf(here('Cartografia/CartografiaOriginal/vegetacion_yunquera.shp'),
                    crs = 25830)


2 Resolución del ejercicio

Se ha dividido el ejercicio en 5 bloques de acuerdo a los objetivos del mismo.

2.1 Flujo de trabajo

Se presenta el flujo de trabajo a continuación. También se puede visualizar directamente aquí (recomendado).


2.2 Mapa de biodiversidad

En esta sección se construirá un mapa donde se representará el Índice de Shanon:

\[ \tag1 H' = - \sum^s_{i=1} p_i \times log_2 p_i \]


donde:

  • i: cada una de las distintas especies

  • s: total de especies distintas

  • \(p_i\): abundancia relativa de cada especie en la comunidad (N individuos de la especie i / N total de individuos).

2.2.1 Descarga de ocurrencias

En primer lugar debemos descargar las ocurrencias de GBIF. En este caso se ha creado un rectángulo en la plataforma web y se han extraído todas las ocurrencias que se encontraban dentro de este rectángulo (correspondiente a la zona de estudio).

A continuación se importan las ocurrencias en R con el nombre de puntos. El siguiente paso es convertir estos puntos en un objeto espacial. Para ello se utiliza la función sf::st_as_sf(). La salida de esta operación (puntos_sf) se transforma al sistema de coordenadas en el que se está trabajando (EPSG: 25830) llamando al objeto ocurrencias. En la Fig. 1 comprobamos que hemos elegido correctamente el rectángulo, y que nuestros puntos se extienden más allá de la zona de estudio asegurando la inclusión de todas las ocurrencias.

puntos <- read.csv(file = here('Cartografia/Tablas/ocurrences.csv'),
                   sep = "\t",
                   quote = "") |> 
  dplyr::select(scientificName,species,decimalLatitude,decimalLongitude) |> 
  drop_na()

puntos_sf <- st_as_sf(puntos,
                      coords = c('decimalLongitude','decimalLatitude'),
                      crs = 4326)

puntos_sf <- st_transform(puntos_sf,
                            crs = 25830)
<center>Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio</center>

Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio


El siguiente paso será delimitar los puntos escogidos al rectánculo mínimo envolvente (bounding box) del área de estudio mediante la función sf::st_crop(). La razón de escoger todos los puntos incluidos dentro del rectángulo y no solamente los que están dentro del área de estudio es que como veremos en la siguiente sección trabajaremos con comunidades ecológicas de forma cuadrangular.

ocurrencias <- st_crop(puntos_sf, yunquera)

Podemos ver el resultado en la Fig. 2, que ascienden a un total de 16906 ocurrencias correspondientes a 1465 especies diferentes.

<center>Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera</center>

Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera


2.2.2 Creación de comunidades ecológicas

Las comunidades ecológicas se asumirán que son los polígonos de la capa de vegetaciones. Para ello, en primer lugar, se crea un identificador único llamado ID y se seleccionan las columnas que más nos interesan. En el siguiente paso se crea una unión (st_join()) entre las ocurrencias y el mapa de vegetación según su intersección. De este modo podemos asignar el ID de la comunidad a cada ocurrencia. Finalmente se transforman las columnas necesarias en data.frame para los posteriores cálculos.

# Crear un ID único
vegetacion <- vegetacion |> 
  cbind(data.frame(ID = sprintf(paste("GID%0",
                                      nchar(nrow(vegetacion)),
                                      "d",
                                      sep=""), 
                                1:nrow(vegetacion)))) |> 
  dplyr::select(ID, D_ARBO_PRE:D_SUEL_PRE, D_USO:COMENTARIO)

# Añadir ID de comunidad ecologica a cada ocurrencia
ocurrencias_id <- ocurrencias |> 
  st_join(vegetacion, join = st_intersects) |> 
  mutate(id_comunidad = factor(ID),
         scientific = scientificName) |> 
  dplyr::select(-ID,-scientificName)

# Crear data frame 
bio <- data.frame(id_comunidad = ocurrencias_id$id_comunidad,
                  scientific = ocurrencias_id$scientific) |> 
  drop_na()
Tabla 1. Encabezado de la tabla de datos de ocurrencias en las distintas comunidades
id_comunidad scientific
GID024 Phoenicurus ochruros (S.G.Gmelin, 1774)
GID031 Abies pinsapo Boiss.
GID043 Erinacea anthyllis Link
GID018 Erinacea anthyllis Link
GID075 Ophrys lutea Cav.
GID100 Abies pinsapo Boiss.


En la Fig. 3 tenemos una representación gráfica de como quedarían estas comunidades con respecto a las ocurrencias, y en la Fig. 4 vemos la densidad de puntos en la zona.

<center>Fig. 3. Comunidades ecológicas en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores</center>

Fig. 3. Comunidades ecológicas en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores


<center>Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera</center>

Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera


2.2.3 Cálculo del índice de Shannon

El último paso de esta sección será calcular el Índice de Shannon (H). Esto lo haremos en cuatro pasos:

  1. Calcular el número de especies por comunidad. En el siguiente cuadro de código se agrupa por comunidad y nombre (group_by), se genera una columna con el número de especies (count()), y se le cambia el nombre a la columna. En la Tabla 2 vemos el encabezado del resultado.
T_num_ind_sp_com <- bio |> 
  group_by(id_comunidad, scientific) |> 
  count() |> 
  rename(num_ind_sp_com = n)
Tabla 2. Encabezado de la tabla del número de especies por comunidad ecológica
id_comunidad scientific num_ind_sp_com
GID001 Avena sativa L. 2
GID001 Thymus mastichina subsp. mastichina 9
GID001 Ulex baeticus subsp. baeticus 5
GID002 Abies pinsapo Boiss. 14
GID002 Asplenium trichomanes L. 1
GID002 Juniperus phoenicea L. 2


  1. Calcular el número de ocurrencias por comunidad. Realizamos un procedimiento similar, pero en este caso solamente agrupamos por comunidad. Vemos el encabezado en la Tabla 3.
T_num_ind_com <- bio |> 
  group_by(id_comunidad) |> 
  count() |> 
  rename(num_ind_com = n)
Tabla 3. Encabezado de la tabla del número de ocurrencias por comunidad ecológica
id_comunidad num_ind_com
GID001 16
GID002 40
GID004 42
GID005 26
GID006 2
GID007 29


  1. El siguiente paso calcula el Índice de Shannon (T_Shannon). Para ello se unen las tablas haciendo un left_join() que incluye todas las columnas de la Tabla 1 de forma que cada comunidad tendrá asociada sus especies, el número de ocurrencias de cada una, y el número total de ocurrencias de la cuadrícula de esa especie. A continuación se crean los coeficientes \(p_i\) y \(\log p_i\) que se utilizarán para el cálculo del índice. En la Tabla 4 vemos la estructura de los datos hasta este punto.
    Finalmente, se calcula el índice de Shannon, se seleccionan solamente las columnas necesarias y se eliminan los duplicados.
T_Shannon <- T_num_ind_sp_com |> 
  left_join(T_num_ind_com,
            by = 'id_comunidad') |> 
  mutate(pi = num_ind_sp_com/num_ind_com,
         logpi = log2(pi)) |> 
  group_by(id_comunidad) |> 
  mutate(shannon = sum(logpi * pi)*-1) |> 
  dplyr::select(id_comunidad, shannon) |> 
  distinct()
Tabla 4. Coeficientes para el cálculo del Índice de Shannon
id_comunidad scientific num_ind_sp_com num_ind_com pi logpi
GID001 Avena sativa L. 2 16 0.1250 -3.000000
GID001 Thymus mastichina subsp. mastichina 9 16 0.5625 -0.830075
GID001 Ulex baeticus subsp. baeticus 5 16 0.3125 -1.678072
GID002 Abies pinsapo Boiss. 14 40 0.3500 -1.514573
GID002 Asplenium trichomanes L. 1 40 0.0250 -5.321928
GID002 Juniperus phoenicea L. 2 40 0.0500 -4.321928


  1. El último paso consiste en dar el valor de Índice de Shannon a cada una de las comunidades. Se añade un paso de reemplazar por 0 en caso de existir NA. Vemos en el webmap la representación final.
# Valor H a cada cuadricula
SF_shannon <- vegetacion |> 
  left_join(y = T_Shannon,
            by = c("ID" = "id_comunidad"))

# Sustituir NA por 0
SF_shannon$shannon <- replace_na(SF_shannon$shannon, 0)


2.3 Segmentación

La segmentación consiste en dividir el territorio en unidades homogéneas en cuanto a unas determinadas características de las variables de entrada. El procedimiento consiste de los siguientes pasos:

  1. Elegir las variables de entrada

  2. Convertir a formato raster y homogeneizar las escalas

  3. Construir un raster virtual

  4. Segmentación del territorio

2.3.1 Variables de entrada

En esta sección se realizarán los pasos 1 y 2 previos. Para ello, construiremos un total de 4 variables (ver cuadro siguiente).

Variables de entrada

Densidad de pinsapo: información sobre la vegetación

Orientaciones: información sobre el estrés abiótico derivado de la situación de solana y umbría

Distancia a canales: información sobre zonas de posible compensación hídrica

Índice de Shannon: información sobre la biodiversidad


Densidad de pinsapo

En primer lugar cargamos la capa, la cual se encuentra en CRS 25830, y la visualizamos con los límites del monte para comprobar que todo es correcto.

SF_pinsapo <- read_sf(here("Cartografia/CartografiaOriginal/inventario_pinsapo.shp"),
                      crs = 25830)


El siguiente paso será crear un raster de densidad. En ese sentido, se utilizará la interpolación inversa a la distancia (IDW).

  1. Primero se transforman los puntos a WGS 1984 ya que las funciones siguientes trabajan este sistema de coordenadas.
SF_pinsapo84 <- st_transform(SF_pinsapo, crs = 4326)
  1. A continuación obtenemos el rectángulo mínimo envolvente(BBox84). Se da una toleracia de 0.005º dado que al reproyectar la imagen se pierde información.
BBox84 <- st_bbox(SF_pinsapo84)


BBox84[c(1,3)] <- BBox84[c(1,3)] + c(-0.005,0)
BBox84[c(2,4)] <- BBox84[c(2,4)] + c(-0.005,0)
  1. Creamos un patrón de puntos con la función spatstat::ppp() utilizando el vector de coordendas X e Y, marks es la variable que utilizaremos para interpolar (Nparc), y en un objeto owin() ponemos las coordendas X e Y del rectángulo mínimo envolvente, que será la extensión de salida.
ppp_pinsapo <- ppp(st_coordinates(SF_pinsapo84)[,1],
                   st_coordinates(SF_pinsapo84)[,2],
                   marks = SF_pinsapo$Nparc,
                   window = owin(BBox84[c(1,3)],
                                 BBox84[c(2,4)]))
  1. Con la función spatstat::idw() realizamos la interpolación y se transforma a raster, que en este caso a diferencia de QGIS no tenemos opción de escoger el número de vecinos de cuáles hacer la interpolación (no obstante, la diferencia del raster de salida es mínima).
R_pinsapo <- ppp_pinsapo |> 
  idw(power = 2, at = 'pixels') |> 
  raster()
  1. Se asigna la proyección y se reproyecta al sistema 25830 cambiando la resolución a 10 metros.
# Asignar proyección
crs(R_pinsapo) <- crs(SF_pinsapo84)

# Transformar a 25830 y asignar resolución espacial
R_pinsapo <- projectRaster(R_pinsapo,
                           crs = 25830,
                           res = 10)
  1. Finalmente se hace una reclasificación de los valores del ráster y se recorta al área de estudio.
# Matriz de entrada para reclasificación
classMatrix <- matrix(c(0,50,0,
                        50,200,100,
                        200,400,300,
                        400,600,500,
                        600,800,700,
                        800,Inf,900),
                      byrow = T,
                      ncol = 3)

# Reclasificación y recorte al área de estudio
R_pinsapo <- R_pinsapo |> 
  reclassify(rcl = classMatrix) |> 
  mask(yunquera)

En el siguiente webmap podemos ver el resultado de los análisis anteriores.


En el resultado vemos que a grosso modo las densidades más altas se encuentran al sureste del monte de Yunquera.


Orientaciones

La segunda de las capas se corresponde con la de orientaciones. En este caso cargamos el modelo digital del terreno (mdt), calculamos las orientaciones en grados con la función raster::terrain(). El siguiente paso será hacer coincidir espacialmente ambos raster con la misma resolución (resample()). Finalmente lo reclasificaremos teniendo en cuenta lo siguiente:

  • 1: zonas de umbría alta (Norte).

  • 2: zonas de umbría (Noreste, Noroeste).

  • 3: zonas medias (Este, Oeste).

  • 4: zonas de solana (Sureste, Suroeste).

  • 5: zonas de solana alta (Sur).

# Cargar capa y cortar a la extensión de Yunquera
mdt <- raster(here("Cartografia/CartografiaOriginal/PNOA_MDT05_ETRS89_HU30_1051_LID.tif"),
              crs = "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs") 

# Cálculo orientaciones
orientaciones <- terrain(mdt,
                         opt = 'aspect',
                         unit = 'degrees') |> 
  resample(R_pinsapo)

# Matriz de clasificación
mat <- matrix(c(-Inf, 22.5,1,
                22.5,67.5,2,
                67.5,112.5,3,
                112.5,157.5,4,
                157.5,202.5,5,
                202.5,247.5,4,
                247.5,292.5,3,
                292.5,337.5,2,
                337.5,Inf,1),
              ncol = 3,
              byrow = TRUE)

# Reclasificación y recorte al área de estudio
R_orientaciones <- orientaciones |> 
  reclassify(rcl = mat) |> 
  mask(yunquera)


Con el webmap y la siguiente tabla a modo informativo vemos que el monte tiene mayor superficie de laderas expuestas a umbría que a solanas.

values(RC_orientaciones) |> 
  tabyl() |> 
  adorn_pct_formatting()
##  values(RC_orientaciones)      n percent valid_percent
##                         1  31553    7.8%         16.7%
##                         2  60604   14.9%         32.0%
##                         3  41307   10.2%         21.8%
##                         4  38322    9.4%         20.2%
##                         5  17696    4.4%          9.3%
##                        NA 216271   53.3%             -


Distancia a canales

El ráster de distancia a canales se ha calculado previamente siguiendo los pasos del guión en QGIS debido a la complicidad de realizar este análisis en R. La capa obtenida se muestra en el siguiente webmap, la cual se encuentra en el sistema de coordenadas correspondiente y con resample() nos aseguramos de que coincida espacialmente con el resto de capas.

R_dist_canales <- raster(here("Cartografia/CartografiaOriginal/Raster_dist_cauces.tif")) |> 
  resample(R_pinsapo) |> 
  mask(yunquera)


Indice de Shannon

El indice de Shannon se ha calculado en la sección 2.2.3. En este caso nos quedaría rasterizar el resultado y remuestrearlo para que coincida espacialmente con el resto de capas.

R_shannon <- SF_shannon |>
  dplyr::select(shannon) |> 
  st_rasterize() |> 
  as("Raster") |> 
  resample(R_pinsapo) |> 
  mask(yunquera)


En el webmap anterior vemos que las zonas de mayor biodiversidad se encuentran en la zona sureste. Se puede interactuar con el mapa poder ver el porcentaje de arbolado y el tipo de formación vegetal en las diferentes áreas.

Guardar las capas

Finalmente se exportan las capas a la carpeta CapasDefRaster.

# writeRaster(R_pinsapo,
#             filename = here("Cartografia/CapasDefRaster/R_pinsapo.tif"),
#             overwrite = T)
# 
# writeRaster(R_orientaciones,
#             filename = here("Cartografia/CapasDefRaster/R_orientaciones.tif"),
#             overwrite = T)
# 
# writeRaster(R_dist_canales,
#             filename = here("Cartografia/CapasDefRaster/R_dist_canales.tif"),
#             overwrite = T)
# 
# writeRaster(R_shannon,
#             filename = here("Cartografia/CapasDefRaster/R_shannon.tif"),
#             overwrite = T)


2.3.2 Raster virtual

El ráster virtual se ha creado en QGIS mediante la herramienta Raster -> Micelanea -> Crear raster virtual.

El resultado se muestra en la Fig. 5.

<center>Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada</center>

Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada


2.3.3 Segmentación del territorio

El último paso de esta sección consiste en realizar una segmentación del monte de Yunquera utilizando el ráster virtual creado en la sección anterior. Para ello, se utilizará la herramienta Orfeo Toolbox en la interfaz de QGIS.

Para llevar esta tarea a cabo, con Orfeo correctamente instalado, se debe ir a Caja de herramientas de Procesos -> OTB -> Segmentation -> Segmentation. Una vez ahí, se seleccionan las opciones de la Fig. 6. La opción “Minimum region size” se probó con varios valores (el resultado debe convertirse a polígono). En el webmap que se verá posteriormente se analizarán las diferencias e influencia de las diferentes capas de entrada en los resultados.

<center>Fig. 6. Parámetros utilizados en la segmentación</center>

Fig. 6. Parámetros utilizados en la segmentación


Para finalizar la tercera tarea de la práctica, se cargan las tres capas creadas con Orfeo en QGIS.

segmentation1500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion1500.shp"))

segmentation2000 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2000.shp"))

segmentation2500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2500.shp"))

En el webmap siguente se puede jugar con las capas para ver el efecto de las distintas segmentaciones en cada una de las capas.


Nota: los bordes de la segmentación no coinciden totalmente con el borde de los ráster. Esto se debe a que las capas se encuentran en crs 25830, pero el webmap se representa en Web Mercator (EPSG: 3857). A modo de análisis no afecta, pero las conclusiones finales se tomarán sobre el Web Map publicado de la siguiente sección.

Los resultados se discuten en las cuestiones finales. No obstante, se hará una breve discusión sobre las tres composiciones de segmentación creadas.

Para analizar los resultados podemos activar en primer lugar Segmentacion (2000), y a continuación Segmentacion (2500) (el orden es importante para su superposición). Tras hacer esto, podemos activar las variables una a una empezando por Densidad pinsapo, Indice de Shannon, Orientaciones y Distancia a canales (m). La segmentación de 2500 píxeles como unidad mínima se tratará como segmentación inicial, y las unidades segmentadas como cantones.

Lo primero que podemos ver es que si segmentamos con unidad mínima de 2000, se crean 11 cantones a mayores. La segmentación inicial parece separar muy bien áreas con alta densidad de pinsapo de áreas con baja densidad. Al aumentar el número de cantones no se ve una clara separación de nuevas áreas con densidad de pinsapo diferentes.

Si activamos el índice de Shannon vemos que la relación no es tan clara como con la variable anterior. Sin embargo, sí que se ve como la segmentación va aproximadamente paralela a zonas de alta biodiversidad para separar zonas de menor biodiversidad. Además, aumentar el número de cantones parece tener una relación algo mayor con esta variable.

Con respecto a la capa de orientaciones, al ser la más irregular, los patrones de segmentación pueden ser más complicados de identificar. Sin embargo, vemos que en zonas donde la densidad de pinsapo y el índice de Shannon (en menor medida) eran aproximadamente uniformes (norte, oeste), la segmentación produce diferencias entre solanas, zonas medias y umbrías bastante claras. El aumento del número de cantones no parece estar muy influenciado por esta variable, aunque al oeste y al norte se pueden ver algunas divisiones más o menos claras.

La variable de distancia a canales también parece tener una gran influencia para separar zonas muy cercanas a canales de zonas muy alejadas. El aumento del número de cantones en la mayoría de los casos se ve muy correlacionado con esta variable.

Si hacemos el mismo procedimiento tomando como capa inicial Segmentación (2000) y la comparamos con Segmentación (1500), vemos que los cantones generados a mayores generan unas irregularidades que a modo de gestión pueden no ser convenientes

Dependiendo del tipo de gestión que llevemos a cabo nos puede interesar más una segmentación u otra. En este monte donde podemos tener preferencias de protección y conservación del pinsapo puede que un tamaño de parcela muy grande no nos interese. En este sentido, el tamaño de píxel mínimo de 2500 crea muchos cantonesde superficie superior a 60 ha por lo que se puede hacer demasiado grande de gestionar para estos objetivos (Fig. 7C). Si activamos en el webmap la capa de segmentación (2000) y densidad pinsapo podemos ver claramente dos cantones con alta densidad de pinsapo situados en la esquina sureste (cantón 69) y otro situado al noroeste de este (cantón 56) que su tamaño puede ser demasiado grande para una gestión más eficiente. Viendo el valor del resto de variables en estos dos cantones, vemos que existe una heterogeneidad notable. Por ello, si activamos segmentación (1500) vemos que estos cantones quedan separados en una mayor número de unidades. En conclusión, se considera que para una gestión forestal del pinsapo nos interesa una segmentación mayor de las zonas de pinsapo que atiendan al resto de variables consideradas. Por ello, para las cuestiones finales se responderá tomando como elección una segmentación con tamaño mínimo de 1500 píxeles.

par(mfrow = c(1,3))
hist(as.vector(st_area(segmentation1500))[as.vector(st_area(segmentation1500))>200]/10000,
     ylab = "",
     xlab = "Superficie (ha)",
     main = "A")

hist(as.vector(st_area(segmentation2000))[as.vector(st_area(segmentation2000))>200]/10000,
     ylab = "",
     xlab = "Superficie (ha)",
     main = "B")

hist(as.vector(st_area(segmentation2500))[as.vector(st_area(segmentation2500))>200]/10000,
     ylab = "",
     xlab = "Superficie (ha)",
     main = "C")
<center>Fig. 7. Histogramas del tamaño de parcela. A) Tamaño mínimo de píxel = 1500; B) tamaño mínimo de píxel = 2000; C) tamaño mínimo de píxel = 2500</center>

Fig. 7. Histogramas del tamaño de parcela. A) Tamaño mínimo de píxel = 1500; B) tamaño mínimo de píxel = 2000; C) tamaño mínimo de píxel = 2500


2.4 Webmap

El cuarto objetivo era crear un webmap y publicarlo. Para ver dicho mapa, hacer click aquí.

2.5 Cuestiones finales

Consideras que el resultado es consistente con la realidad que conoces del terreno? Justifica tu respuesta.

Teniendo en cuenta las variables consideradas sí que es consistente el resultado. No obstante, si superponemos la segmentación junto a una imagen aérea podemos ver que al no tener en cuenta los caminos o pistas, la segmentación ignora estos elementos del paisaje que son clave en la creación de cantones para una gestión más sencilla y eficiente. Por ello, a efectos prácticos puede perder consistencia. Por otro lado, si nos fijamos en cuando a homogeneidad de las masas forestales en cuanto a su densidad (en la imagen aérea) podemos ver una alta correlación.

Cuál de las capas (variables) utilizada aparenta tener más influencia en el resultado final? Por qué crees que pasa esto?

Densidad pinsapo > distancia a canales > índice de Shannon > orientación.

En las zonas de cambio de densidad de pinsapo, el resto de variables parecen ser también más o menos homogéneas, y en donde cambia la densidad de pinsapo, es donde se producen cambios mayores en la densidad de los píxeles dado que esta variable produce unos cambios de densidad mayores al resto.

Para qué aspectos de la gestión del territorio crees que sería de utilidad la cartografía obtenida?

Principalmente para una gestión forestal del pinsapo ya que tenemos bien separados los cantones de mayor densidad que los de menor densidad. Si añadimos un estudio de fauna, podríamos analizar qué especies están presentes en el monte y por qué unas zonas tienen una biodiversidad mucho mayor que otras para, en medida de lo posible y sin perjudicar otros usos que pueda tener el monte que no se quieran perder, aumentar la biodiversidad en otras zonas.

Cómo mejorarías el resultado? Incorpora alguna sugerencia sobre mejoras técnicas y/o nuevas variables a incorporar.

Para añadir nuevas variables y mejorar el resultado es indispensable saber qué tipo de gestión de llevará a cabo. En este caso se considera la gestión forestal del pinsapo:

Para realizar gestión forestal, es interesante aprovechar la red viaria disponible en el monte ya que se definen límites más claros para el gestor forestal y porque son aprovechados para la extracción de madera. Además, si cortamos el MDT a la extensión del monte vemos que existe un gradiente altitudinal de 1000 metros, que va a influir mucho en el desarrollo de la vegetación. De este modo, añadir una variable de distancia a caminos y la altitud podría ser una buena idea para una gestión forestal. Finalmente podríamos suavizar los bordes de los cantones.

Información de la sesión

## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3     spatstat_3.0-2         spatstat.linnet_3.0-3 
##  [4] spatstat.model_3.0-2   rpart_4.1.19           spatstat.explore_3.0-5
##  [7] nlme_3.1-160           spatstat.random_3.0-1  spatstat.geom_3.0-3   
## [10] spatstat.data_3.0-0    raster_3.6-3           sp_1.5-1              
## [13] stars_0.5-6            abind_1.4-5            here_1.0.1            
## [16] terra_1.6-17           rgbif_3.7.3            knitrdrawio_0.2.2     
## [19] knitr_1.40             sf_1.0-9               forcats_0.5.2         
## [22] stringr_1.4.1          dplyr_1.0.10           purrr_0.3.5           
## [25] readr_2.1.3            tidyr_1.2.1            tibble_3.1.8          
## [28] ggplot2_3.4.0          tidyverse_1.3.2        janitor_2.1.0         
## [31] leafpop_0.1.0          mapview_2.11.0        
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.1            uuid_1.1-0              backports_1.4.1        
##   [4] systemfonts_1.0.4       lwgeom_0.2-9            plyr_1.8.8             
##   [7] lazyeval_0.2.2          splines_4.2.1           crosstalk_1.2.0        
##  [10] leaflet_2.1.1           digest_0.6.30           htmltools_0.5.3        
##  [13] leaflet.providers_1.9.0 fansi_1.0.3             magrittr_2.0.3         
##  [16] tensor_1.5              googlesheets4_1.0.1     tzdb_0.3.0             
##  [19] modelr_0.1.10           svglite_2.1.0           timechange_0.1.1       
##  [22] spatstat.sparse_3.0-0   colorspace_2.0-3        rvest_1.0.3            
##  [25] rgdal_1.6-2             haven_2.5.1             xfun_0.34              
##  [28] leafem_0.2.0            crayon_1.5.2            jsonlite_1.8.3         
##  [31] brew_1.0-8              glue_1.6.2              polyclip_1.10-4        
##  [34] gtable_0.3.1            gargle_1.2.1            webshot_0.5.4          
##  [37] scales_1.2.1            oai_0.4.0               DBI_1.1.3              
##  [40] Rcpp_1.0.9              viridisLite_0.4.1       units_0.8-0            
##  [43] klippy_0.0.0.9500       proxy_0.4-27            stats4_4.2.1           
##  [46] htmlwidgets_1.5.4       httr_1.4.4              ellipsis_0.3.2         
##  [49] pkgconfig_2.0.3         farver_2.1.1            sass_0.4.2.9000        
##  [52] dbplyr_2.2.1            deldir_1.0-6            utf8_1.2.2             
##  [55] conditionz_0.1.0        tidyselect_1.2.0        labeling_0.4.2         
##  [58] rlang_1.0.6             munsell_0.5.0           cellranger_1.1.0       
##  [61] tools_4.2.1             cachem_1.0.6            cli_3.4.1              
##  [64] generics_0.1.3          broom_1.0.1             evaluate_0.19          
##  [67] fastmap_1.1.0           yaml_2.3.6              goftest_1.2-3          
##  [70] fs_1.5.2                satellite_1.0.4         whisker_0.4            
##  [73] xml2_1.3.3              compiler_4.2.1          rstudioapi_0.14        
##  [76] png_0.1-7               e1071_1.7-12            spatstat.utils_3.0-1   
##  [79] reprex_2.0.2            bslib_0.4.1             stringi_1.7.8          
##  [82] highr_0.10              lattice_0.20-45         Matrix_1.5-1           
##  [85] classInt_0.4-8          vctrs_0.5.0             pillar_1.8.1           
##  [88] lifecycle_1.0.3         jquerylib_0.1.4         data.table_1.14.4      
##  [91] R6_2.5.1                KernSmooth_2.23-20      codetools_0.2-18       
##  [94] MASS_7.3-58.1           assertthat_0.2.1        rprojroot_2.0.3        
##  [97] withr_2.5.0             mgcv_1.8-41             parallel_4.2.1         
## [100] hms_1.1.2               grid_4.2.1              prettydoc_0.4.1        
## [103] class_7.3-20            rmarkdown_2.18          snakecase_0.11.0       
## [106] googledrive_2.0.0       lubridate_1.9.0         base64enc_0.1-3